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Abstract The quantum XY model shows a Berezinskii-Kosterlitz-Thouless (BKT) tran- 
sition between a phase with quasi long-range order and a disordered one, like 
the corresponding classical model. The effect of the quantum fluctuations is to 
weaken the transition and eventually to destroy it. However, in this respect the 
mechanism of disappearance of the transition is not yet clear. In this work we 
address the problem of the quenching of the BKT in the quantum XY model 
in the region of small temperature and high quantum coupling. In particular, 
we study the phase diagram of a 2D Josephson junction array, that is one of 
the best experimental realizations of a quantum XY model. A genuine BKT 
transition is found up to a threshold value g* of the quantum coupling, beyond 
which no phase coherence is established. Slightly below g* the phase stiffness 
shows a reentrant behavior at lowest temperatures, driven by strong nonlinear 
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quantum fluctuations. Such a reentrance is removed if the dissipation effect of 
shunt resistors is included. 



1. 



Introduction 



Two-dimensional (2D) arrays of Josephson junction (JJA) are one of the 
best experimental realizations of a model belonging to the XY universality 
class and permit to check and study a variety of phenomena related to both the 
thermodynamics and the dynamics of vortices. In these systems a Berezinskii- 
Kosterlitz-Thouless (BKT) transition [1] separates low-temperature supercon- 
ducting (SC) state from the normal (N) state, the latter displaying no phase co- 
herence [2]. At nanoscale size of the junctions, the quantum fluctuations of the 
superconducting phases cause new interesting features. These appear to be the 
consequence of the non-negligible energy cost of charge transfer between SC 
islands. Indeed small capacitances are involved and the phase and charge are 
canonically conjugate variables. A relevant effect is the progressive reduction 
of the SC-N transition temperature, an example of which is shown in Fig. 1.1, 
where experimental data [3] are compared with semiclassical results [4]. 
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Figure 1.1. Experimental phase diagram for a JJA on triangular (left panel) and square (right 
panel) lattice [2]. The lower temperatures refer to the superconducting phase. The lines report 
the semiclassical result of the PQSCHA [4], in the limit R s > Rq = h/(2e) 2 with r) = 1CT 2 
(solid), and for strong dissipation, Rq/Rs = 3 (dashed). 



Recently, fabricated arrays of nanosized junctions, both unshunted [3] and 
shunted [5], have given the opportunity to experimentally approach the quan- 
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turn (zero-temperature) phase transition. However, the mechanism of suppres- 
sion of the BKT in the neighborhood of the quantum critical point and its 
connection with the observed reentrance of the array resistance as function of 
the temperature is not yet clear [2, 3, 6]. In this paper we study the SC-N phase 
diagram by means of path-integral Monte Carlo (PIMC) [7, 8] simulations fo- 
cusing the attention on the region of strong quantum fluctuations, in order to 
investigate their role in suppressing the BKT transition. 

2. The model 

The Josephson junction two-dimensional array (JJA) on the square lattice is 
modeled by a quantum XY model with the following action: 



where tpij = <p% — <fj is the phase difference between the Josephson phases 
on the ith and the jth neighboring superconducting islands. The capacitance 
matrix reads 



where Co = f] C and C are, respectively, the self- and mutual capacitances 
of the islands, and d runs over the vector displacements of the z = 4 nearest- 
neighbors. The standard samples of JJA are well decribed by the limits rj <C 1, 
while for the granular films the opposite limits rj 3> 1 is more appropriate. 
The quantum dynamics of the system is determined by the Coulomb interac- 
tion between the Cooper pairs. This is described by the kinetic term through 
Josephson relation ipi = 2eVi. The Josephson coupling is represented by the 
cosine term, the latter favoring Cooper-pair tunneling across the junctions. 

The quantum fluctuations are ruled by the quantum coupling parameter 
g = \/E c / Ej , where E c = (2e) 2 /2C is the characteristic charging en- 
ergy (for t] <C 1). In the following we use the dimensionless temperature 
t = k-gT/Ej. In our model Eq. (1.1) we assume the presence of very weak 
Ohmic dissipation due to the currents flowing to the substrate or through shunt 
resistances [5], which reflects into the prescription to consider the phase as an 
extended variable [2]. Apart from this, dissipative effects are negligible pro- 
vided that the shunt resistance R s 3> R Q g 2 /(2nt), where R Q = h/(2e) 2 
is the quantum resistance; for smaller R s an explicit dissipative contribution 
should be added to the action (1.1), e.g. in the form of the Caldeira-Leggett 




(1.1) 




(1.2) 



term [2, 4], 




o 



o 



(1-3) 



4 



1.15 : 

1.1 : 

1.05 : 

1 - 

0.95 : 

0.9 : 

0.85 : 





<r. 




>o 




41 














II • 


II 


II 


II 


II 


II 


6* 


6* 


60 ! 


*« 







0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 

t 



Figure 1.2. Fitting parameter A(t) for different values of g ranging from 0.5 to 3.0. The BKT 
transition temperature is obtained by the condition A(t BKT ) = 1. 



resulting in a decrease of quantum fluctuations. The two situations are the cases 
of the experiments in Ref. [3] and Ref. [5], respectively, where an increasing of 
the BKT transition temperature was found for increasing dissipation. This can 
be easily understood taking into account that the dissipative term (1.3) results 
from the coupling of the phase ipi with environmental variables (the degrees of 
freedom of the dissipative bath), constituting an implicit measurement of ipi . 

3. Numerical simulations 

The numerical data for the BKT transition temperature are obtained using 
PIMC simulations onLxI lattices (up to L = 96) with periodic boundary 
conditions. Thermodynamic averages are obtained by MC sampling of the par- 
tition function after discretization of the Euclidean time u G [0, j3H\ in P slices 
h(3/P , where P is the Trotter number. However, the actual sampling is made 
on imaginary-time Fourier transformed variables on a lattice using the algo- 
rithm developed in Ref. [7]. The move amplitudes are independently chosen 
and dynamically adjusted for each Fourier component; this procedure turns out 
to be very efficient to reproduce the strong quantum fluctuations of the paths 
in the region of high quantum coupling g. Indeed, test simulations with the 
standard PIMC algorithm showed serious problems of ergodicity, though even- 
tually giving the same results. The approach through Fourier PIMC becomes 
more and more suitable and effective when dissipation is inserted, because the 
non-local action (1.3) becomes local in Fourier space. Additional details about 
the numerical method are given in the Appendix. An over-relaxation algo- 
rithm [9] over the zero-frequency mode has also been implemented in order to 
effectively reduce the autocorrelation times. 
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A very sensitive method to determine the transition temperature is provided 
by the scaling law of the helicity modulus T (a quantity proportional to the 
phase stiffness), 



Ej V dk 2 Q , kQ=Q 

which measures the response of the free energy F (per unit volume) when a 
uniform twist ko along a fixed direction u is applied to the boundary conditions 
(i.e., <fi — > (fii + kou-i, with the unitary vector u). 

The PIMC estimator for T is easily obtained, in analogy to that of Ref. [10], 
by derivation of the path-integral expression of the partition function (see Ap- 
pendix). Kosterlitz's renormalization group equations provide the critical scal- 
ing law for the finite-size helicity modulus T^: 

TL(t BKT ) = 2 / 1 \ 

t BKT vr V + 2HL/L Q )) ' 1 

where L$ is a non-universal constant. Following Ref. [11], the critical temper- 
ature can be found by fitting Ti(t)/t vs L for several temperatures according 
to Eq. (1.5) with a further multiplicative fitting parameter A{t). In this way, 
the critical point can be determined by searching the temperature such that 
i(t BKT ) = 1, as illustrated in Figs 1.2 and 1.3. Using this procedure the crit- 
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Figure 1.3. Size scaling of the helicity modulus Tl at the transition temperature. Symbols 
are PIMC data and the dashed-lines are the one-parameter fit with Eq. (1.5), i.e. with A(t) = 1. 
Left panel: g = and t = 0.892 [L = 0.456(6)]; right panel: g = 3.4 and t = 0.25 
[Lo = 3.32(3)]. The insets show A(t) for different temperatures, using the two-parameter fit 
(see text). 
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Figure 1.4. Helicity modulus Tl divided by the best fit with the expression (1.5) for g = 3.4 
and different temperatures: A, 0> an d D correspond to t = 0.2, 0.25, 0.3, respectively. 

ical temperature can be determined with excellent precision. For instance in 
the classical case we get t BKT (g=0) = 0.892(2), in very good agreement with 
the most accurate results from classical simulations [12]. Also in the regime 
of strong quantum coupling, g = 3.4, the PIMC data for T^(t BKT =0.25) 
are very well fitted by Eq. (1.5), as shown in Fig. 1.4. Moreover, this fig- 
ure points out the sensitivity of this method to identify t BKT : at temperature 
higher (lower) than the critical one the helicity modulus decreases (increases) 
much faster with L than Ti(t BKT ). At higher values of the quantum coupling, 
g > g*, the helicity modulus scales to zero with L — > oo and P — > oo at any 
temperature [8]. 

4. Results 

We have found significant differences with the standard BKT theory. In 
the regime of strong quantum fluctuations, a range of coupling values, 3.2 < 
g < g*, is found in which the helicity modulus displays a non-monotonic 
temperature behavior. 

In Fig. 1.5, Ti(t) is plotted for different values of g on the 8x8 cluster: up 
to g = 3.0 it shows a monotonic behavior similar to the classical case, where 
thermal fluctuations drive the suppression of the phase stiffness. In contrast, for 
g = 3.4, the helicity modulus is suppressed at low temperature, then it increases 
up to t ~ 0.2; for further increasing temperature it recovers the classical- 
like behavior and a standard BKT transition can still be located at t ~ 0.25 
(Figs. 1.3 and 1.4). A reentrance of the phase stiffness was already found for 
a related model in Ref. [10], but the authors concluded that the drop of the 
helicity modulus at lowest temperatures was probably due to the nniteness of 
the Trotter number P . 
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Figure 1.5. Temperature behavior of the helicity modulus Tg (t) on a 8 x 8 lattice, for different 
values of g. The data are results from the Trotter extrapolation. 

Systematic extrapolations in the Trotter number and in the lattice size have 
been done, and presented in Figs. 1.6 and 1.7 for g = 3.4 , in order to ascertain 
this point. In particular, we did not find any anomaly in the finite-P behavior: 
the extrapolations in the Trotter number appear to be well-behaved, in the ex- 
pected asymptotic regime 0(1/ P 2 ) [13], for P > 60 (Fig. 1.6). Moreover, 
the extrapolation to infinite lattice-size shown in Fig. 1.7 clearly indicates that 
Tl scales to zero at t = 0.1, while it remains finite and sizeable at t = 0.2. 
therefore, the outcome of our analysis is opposite to that of Ref. [10], i.e., we 
conclude that the reentrant behavior of the helicity modulus appears to be a 
genuine effect present in the model, rather than a finite-Trotter or finite-size 
artifact. 

In order to understand the physical reasons of the reentrance observed in the 
phase stiffness, we have studied the following two quantities: 

A = ( cos tpij(u) ) , (1.6) 
Al = ((iPijiu)-^) 2 ) , (1.7) 

with 

1 fW 

and ij nearest-neighbor sites. The first quantity A , is a measure of the total 
(thermal plus quantum) short-range fluctuations of the Josephson phase and 
its maximum occurs where the overall fluctuations are lowest. The second 
quantity represents instead the pure-quantum spread of the phase difference 
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Figure 1.6. Trotter-number extrapolation of Tl for g — 3.4. Two series of data for t = 0.1 
(•) and 0.2 (■) are reported, for four different lattice sizes: from the top to the bottom L = 
8, 10, 12, 14. The lines are weighted quadratic fits. 
0.3 



0.2 - 



0.1 







■ ' t = 0.2 












/ 

/ 

f^t = 0.1 ~ 



0.04 



0.08 

1/L 



0.12 



0.16 



Figure 1. 7. Finite-size scaling of the helicity modulus Tl for g = 3.4 at fixed P = 101. The 
lines are guides for the eye. 



between two neighboring islands and has been recently studied in the single 
junction problem [14]; more precisely, A^, measures the fluctuations around 
the static value (i.e., the zero-frequency component of the Euclidean path), it 
is maximum at t = and tends to zero in the classical limit, i.e., (g/t) —> 0. 

The quantities (1.6) and (1.7) on a 8 x 8 lattice are compared in Fig. 1.8 
for two values of the quantum coupling, in the semiclassical (g = 1.0) and in 
the extreme quantum (g = 3.4) regime. In the first case A = (cosipij(u)) 
decreases monotonically by increasing t and the pure-quantum phase spread 
A^, shows a semiclassical linear behavior which is correctly described by the 
PQSCHA. At variance with this, at g = 3.4, where the reentrance of Y(i) 
is observed, A = (cos (p%j(u)) shows a pronounced maximum at finite tem- 
perature. Besides the qualitative agreement with the mean-field prediction of 
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Ref. [15], we find a much stronger enhancement of the maximum above the 
t = value. This remarkable finite-i effect can be explained by looking at 
A^,(<7=3.4) (Fig. 1.8, top-right panel). In fact its value is an order of mag- 
nitude higher than the one in the semiclassical approximation and, notably, 
it is strongly suppressed by temperature in a qualitatively different way from 




Figure 1.9. Qualitative picture of the behavior of the total mean square fluctuation of the 
nearest-neighbor phase difference. The straight line refer to the classical fluctuations. The lower 
(dashed) curve refers to the low-coupling case. The upper curve represents strongest couplings. 
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A^(<7=1.0): i.e., the pure-quantum contribution to the phase fluctuations mea- 
sured by Az, decreases much faster than the linearly rising classical (thermal) 
one resulting in a global minimum of the total fluctuations. The qualitative 
behavior of the total mean square fluctuation of the nearest-neighbor phase 
difference vs. temperature is sketched in Fig. 1.9. This single-junction effect, 
in a definite interval of the quantum coupling (3.2 < g < 3.4), is so effec- 
tive to drive the reentrance of the phase stiffness. As for the transition in the 
region of high quantum fluctuations and low temperature, the open symbols 
in Figs. 1.10 and 1.11 represent the approximate location of the points (£, g) 
where T(i) becomes zero within the error bars: in their neighborhood we did 
not find any BKT-like scaling law. This fact opens two possible interpretations: 
(i) the transition does not belong to the XY universality class; (it) it does, and 
in this case the control parameter is not the (renormalized) temperature, but a 
more involute function of both t and g. 

5. The phase diagram 

Let us now describe our resulting phase diagram [8], displayed in Fig. 1.10 
together with the semiclassical results valid at low coupling [4] and in Fig. 1.11 
together with the experimental data [3] and our first PIMC outcomes including 
the dissipation effect. 
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Figure 1.10. Phase diagram for a square lattice Josephson junction arrays in the limit Rs 
R Q = h/(2ef with rj = 1CT 2 . The symbols are our PIMC results and the line reports the 
semiclassical result of the PQSCHA [4]. 



At high temperature, the system is in the normal state with vanishing phase 
stiffness and exponentially decaying phase correlations, (ip\(p- s ). By lowering 
t, for g < g* ~ 3.4, the system undergoes a BKT phase transition at t BKT (g) 
to a superconducting state with finite stiffness and power-law decaying phase 
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correlations. When g is small enough (semiclassical regime), the critical tem- 
perature smoothly decreases by increasing g and it is in remarkable agreement 
with the predictions of the pure-quantum self-consistent harmonic approxima- 
tion (PQSCHA) [4]. For larger g (but still g < g*) the semiclassical treatment 
becomes less accurate and the curve t BKT (g) shows a steeper reduction, but 
the SC-N transition still obeys the standard BKT scaling behavior. Finally, 
for g > g* a strong quantum coupling regime with no sign of a SC phase is 
found. Surprisingly, the BKT critical temperature does not scale down to zero 
by increasing g (i.e., t BKT (g*) ^ 0): by reducing the temperature in the region 
3-2 < g < g* ', phase coherence is first established, as a result of the quenching 
of thermal fluctuations, and then destroyed again due to a dramatic enhance- 
ment of quantum fluctuations near t = 0. This is evidenced by a reentrant 
behavior of the stiffness of the system, which vanishes at low and high t and it 
is finite at intermediate temperatures. The open symbols in Fig. 1.10 mark the 
transition between the finite and zero stiffness region when t is lowered. 

When the interaction with a heat bath given by Eq. (1.3) is present through 
a variable shunt resistance, the quantum phase fluctuations are decreased by 
the dissipation so that the BKT transition temperature rises. This is well repro- 
duced by PQSCHA approach [4] and turns out to be in qualitative agreement 




Figure 1,11. Phase diagram for a square lattice Josephson junction arrays. Crossed circles are 
the experimental data [2]. The other filled symbols are our PIMC results, in the limit _Rs 3> 
Rq — h/(2e) 2 (diamonds), for Rq/Rs = 0.25 (squares) and Rq/Rs = 0.5 (circles), all with 
r\ — 10 -2 . The solid line reports the semiclassical non dissipative result of the PQSCHA [4]. 

with recent experiments [5]. As shown in Fig. 1.11, the presence of interme- 
diate (Rq/Rs = 0.25) and high (Rq/Rs = 0.5) dissipation causes the dis- 
appearance of the reentrance, in agreement with the experiments of Takahide 
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et al. [5]. Their data do not show any anomalous reentrant behavior like those 
ones of the experiments on unshunted samples in Ref. [3]. Further investiga- 
tions are needed to answer the question about the nature of the transition in 
the reentrance zone and to determine at which values of the dissipation the 
reentrant behavior is washed out. 



6. Summary 

In summary, we have studied a model for a JJA in the quantum-fluctuation 
dominated regime by means of Fourier path-integral Monte Carlo simulations. 
The BKT phase transition has been followed increasing the quantum coupling 
g up to a critical value g* ~ 3.4 where t BKT ~ 0.25; above g* no traces of BKT 
critical behavior have been observed. Remarkably, in the regime of strong 
quantum coupling (3.2 < g < 3.4) phase coherence is established only in a 
finite range of temperatures, disappearing at higher T, with a genuine BKT 
transition to the normal state, and at lower T, due to a nonlinear quantum 
mechanism. This effect is destroyed by the presence of a relevant dissipation. 
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Appendix: PIMC in the Fourier space for JJA 

In this appendix we derive the discretized expressions of the observables to be measured by 
means of the simulations based on the Metropolis algorithm. The first step is to discretize the 
Euclidean time in P slices, m = (h/3/P)l, and to perform a discrete Fourier transform 

= Jl>M <si * ,a , (1.A.1) 
i=i 

where tpij = ipi{u{), and (p i>a satisfies the conditions of periodicity, i.e. ip i>a — ipi >a +p, 
and reality of iptj, i.e. <^i iQ = Thus the discretized path can be written as a sum over 

different frequency components 



N r 

2nla , . 2nla 
Oi, a cos — — I- bi, a sin 



p 

(1.A.2) 

where (pi is the zero-frequency component of the Euclidean path Eq. (1.8), and we choose an 
odd Trotter number P — 2N + 1. The advantage of using expression (1.A.2) is twofold. On 
the one hand the dissipative term of the action (1.3) that is non-local in time becomes diagonal. 
On the other hand the sampling can be performed independently on each frequency component, 
dynamically adjusting each move amplitude. This method make the sampling very effective 
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especially in the region of strong quantum fluctuations, where the main contribution to the path 
comes from {cii iCt } and {bi, Q }. 

Using expression (1.A.2), the JJA action (1.1) plus the dissipative term (1.3) reads 



t , s OTP 

S l<P] = Z Z T ijA a i,aaj, a + h,ofij, a ) + -p 1 ^ Z ( X ~~ COSt Pij,l) , (1-A.3) 



ij a — 1 (ij) 1 = 1 

where (fij,i = <Pi,l—<Pj,l, and the "kinetic" matrix is 



T- ■ 



(1.A.4) 



and Kij ta is the discrete FT of the dissipative kernel matrix Kij(u—u'). All the macroscopic 
thermodynamic quantities are obtained through the estimators generated from the discretized 
action (1.A.3). For instance the estimator of the helicity modulus T is obtained applying the 
definition (1.4) to the discretized free energy per unit volume 
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with the normalization constant C. Eventually we get 
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